MULTI-CONFIGURATIONAL SYMMETRIZED PROJECTOR 
QUANTUM MONTE CARLO METHOD FOR LOW-LYING EXCITED 
STATES OF THE HUBBARD MODEL Q 



Bhargavi Srinivasan 3 , S. Ramasesha 3 ' 5 and H. R. Krishnamurthy 4,5 f] 

3 Solid State and Structural Chemistry Unit 

^Department of Physics 
Indian Institute of Science, Bangalore 560 012, India 

5 Jawaharlal Nehru Centre for Advanced Scientific Research 
Jakkur Campus, Bangalore 560 064, India 

PACS numbers: 71.10Fd, 71.20Ad, 71.20Hk 



Contribution no. 1179 from the Solid State and Structural Chemistry Unit 
2 e-mail: bhargavi@sscu.iisc. ernet. in, ramasesh@sscu.iisc. ernet. in, 
hrkrish@physics . iisc . ernet . in 



1 



ABSTRACT 



We develop a novel multi-configurational Symmetrized Projector Quantum Monte Carlo 
(MSPQMC) method to calculate excited state properties of Hubbard models. We compare 
the MSPQMC results for finite Hubbard chains with exact results (where available) or 
with Density Matrix Renormalization Group (DMRG) results for longer chains. The 
energies and correlation functions of the excited states are in good agreement with those 
obtained from the alternate techniques. 

1 Introduction 

In recent years, numerical many-body techniques which have proved reliable for the study 
of the Hubbard model are the Projector Quantum Monte Carlo (PQMC) method]!], 0] 
and the Density Matrix Renormalization Group (DMRG)[[| method (for 1-D and quasi 

1- D systems). These methods have been predominantly ground state techniques. The 
usual procedure for obtaining excited states of a many-body Hamiltonian in exact diag- 
onalization schemes is to exploit the symmetries of the system and block-diagonalize the 
Hamiltonian [|]]. The lowest few eigenvalues in each block can then be computed using 
standard numerical procedures. However, exact diagonalization schemes are limited to 
rather small system sizes, unlike either the DMRG scheme or the PQMC scheme. Im- 
plementation of symmetries in the latter schemes turns out to be nontrivial although 
desirable. 

While the DMRG method could yield a few low-lying states, low-lying excited states 
of a chosen symmetry were inaccessible from this technique, until its recent extension 
to incorporate crucial symmetries of a system^], which enables the method to target 
excited states as low-lying states in subspaces of a given irreducible representation of the 
symmetry group of the given system. 

In contrast, the PQMC method has exclusively been a ground state technique for 
fermionic systems. Hitherto, it was not feasible to obtain even the ground state of the 

2- D Hubbard model for arbitrary filling because of the open-shell structure of the non- 
interacting ground state ||. Here, we present a novel multi-configurational symmetrized 
PQMC (MSPQMC) technique which makes it possible, for the first time, to obtain ener- 
gies of excited states of the Hubbard Hamiltonian. The technique is also applicable to the 
ground state of open shell systems. We use the recently developed symmetrized PQMC 
method^ to improve the Monte Carlo estimates of properties in the targetted state. In 
the next section, the formulation of the MSPQMC method will be presented. In section 3, 
we demonstrate the method by applying it to Hubbard chains. We also discuss numerical 
issues associated with the MSPQMC procedure. 
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2 The MSPQMC Method 



The single band Hubbard Hamiltonian H for a system of N sites, may be written as||, 

N 

H = H + Hi = ~(Y1 %«L%<t + h.c.) + U^2h^hi h (1) 

(ij),cr 1=1 

where the symbols have their usual meanings. 

Using the projection ansatz, the lowest eigenstate, >, in a given irreducible sym- 
metry subspace T, of H, can be projected from a trial wavefunction |0 r > as 



e-P H \6 T > 



|V>o >= lim 

ft-* 00 sl< r |e" 2 ^|0 r > 

provided \<p v > has a nonzero projection on to |-0q >. The trial wavefunction \cf) T > is 
usually formed from the molecular orbitals (MO) obtained as eigenf unctions of the non- 
interacting part, Hq, °f the full Hamiltonian. When the non- interacting ground state 
of a given system is a closed-shell state, the trial wavefunction \<p v > for obtaining the 
interacting ground state can be chosen as a single nondegenerate electronic configuration 
in the MO basis. However, to obtain an excited state as the lowest eigenstate of the 
interacting model within a desired symmetry subspace from the projection ansatz, it is 
usually necessary to choose |0 r > to be a specific linear combination of degenerate excited 
MO-configurations as the trial wavefunction. Such a linear combination can be obtained 
by operating with the group theoretic projection operator[|l(J for the desired irreducible 



representation on a single excited MO-configuration. In order to fix the total spin of the 
target state, we use the L6wdin|| projection operator to project out the desired spin 
state from the trial configuration transforming as T. The projection procedure in eqn. 
(0) conserves the symmetry of the initial state and hence projects out the lowest energy 
excited state of the interacting model with that symmetry from the trial state. The trial 
state |0 r > in general takes the form, 

l^>=E C M> i l# >= IC > l#,-«r >. ( 3 ) 
3=1 

where p is the number of degenerate MO-configurations in the symmetry adapted starting 
wavefunction. An MO-configuration with M a fermions of spin a in second quantized form 
can be written as 

M CT N 

ic>= n(£(*£Wk)io> (4) 

m=l i=l 
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where 3>£ r is an N x M a sub- matrix of the MO coefficients whose row index, i, labels sites 
and the column index, m, labels the MOs occupied by electrons of spin a, in the MO- 
configuration labelled by j and T. The overlap of any two MO-configurations expressed 
in this form is given by 



< CrWv >= det 



(5) 



In the PQMC method for the Hubbard model, the projection operator exp(—/3H) is 
Trotter decomposed as (exp(—ArH)) L with L imaginary time slices of width At (/3 = 
L x At). This is followed by a discrete Hubbard-Stratanovich (H-S) transformation ]lj ] 



of the on-site interaction Hamiltonian at each site i and each time-slice /, in terms of 
Ising-like fields, sn, leading to the result 



e 



-AtH 



^X^s^X^il^i) (6) 



fa} 



— At - AtU —At * 
X a (l, Sj) = exp[^—H ] 22 exp[( a \ s a n ia — }exp[^— H ] (7) 

where the summation is over all possible iV-vectors S/ whose i th components correspond 
to the H-S field, s#, ( a is +1 (-1) for electrons with | (J.) spin and the H-S parameter 



A = 2arctanh^tanh(ATU /A). Thus, 

e"^ = E^(W)WMW) = Y,W({s}) (8) 

W^(M) = X a (L,s L )...X a (l, Sl ) (9) 

The action of each of the terms in the summation in eqn. (^j) on a trial state of the form 
in eqn. (|j) can be obtained as the left multiplication of the iV x M G matrix by an 
N x N matrix, B CT (/, S;), given by 

B a (l, = b b lCT (/, s ; )b (10) 

The matrix b is given by exp[— K], with Ky = — 4r-ty ■ The matrix b lcr (/, S/) is diagonal 
with elements 5y|exp[Co-Asa — ^]. 

To obtain the expectation value of an operator O in the targetted state, when the 
trial state is a single MO-configuration, < O >= (< ?/> r |0|'?/> r >)/(< ip r \ip r >) we define 
states \R(1,{sr}) > and < L(l,{s L })\. The former is obtained by projecting the trial 
wavefunction through the right Ising lattice {sr} formed by time-slices 1 through I, while 
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the latter is obtained by projecting its transpose through the left Ising lattice {sl} time- 
slices L through L — I + 1, 

\^>^j2\R r (i,{ SR })> = E^(W))l0 r > (ii) 

{ s r} { s r} 

< ^| « E < ^ W})| = E < (12) 

This allows us to express < O > as a weighted average, X) a; r ({s})O r ({s}), with weights 

{4 

given by, 

= 

E<0 r l^(W)l0 r > 

{4 

n rvr ,x < £ r (Z,{s L })|0|i? r (Z,W})> 

U >J <L r (U*4)l^ r (U*4)> ' 1 J 

For example, if the operator O is the single-particle operator al a a>i a , O r ({s}) takes the 
form, 

O r ({s}) = ((R^(l,{s}))(L^(l,{s})R^(l,{s})- 1 (L^(l,{s}))) ki . (15) 

which is the kl th element of the single-particle Green function, G^(l, If we weight 

average the property over all the Ising-configurations, we would obtain the expectation 
value of that property in the targetted state, exact to within Trotter error. However, 
exhausting all Ising-configurations in an averaging procedure is impractical and the de- 
nominator in the eqn. ([13]) cannot be known explicitly. Therefore, we resort to an 
importance sampling Monte Carlo (MC) estimation in which a knowledge of the ratio 
of weights, r, for any two configurations {s'} and {s}, uj r ({s'})/u) r ({s}) is sufficient for 
obtaining property estimates. The ratio, r is given by the ratio of inner products of the 
right and left projected states (using eqn. (R)) for the two Ising-configurations, 



det(Ll(l,{s' L })RUl,W R })) 
K l det{m,{s L })K(l,{s R })) 



Ising-configurations are generated by sequential single spin-flips through the lattice, exam- 
ining each site at a given time slice, I, before proceeding to the next. This allows efficient 
computation of the ratio, r. Using the heat bath algorithm, the new configuration is 
accepted or rejected with a probability r/(l + r). 

The above procedure can be extended to a multi-configuration trial function (p > 1 
in eqn.(PD) and the ratio of weights for Ising-configurations {s'} and {s} takes the form, 

lu v ({s'}) <L r (l,{s'})\R r (l,{s'})> 



^({s}) <L^l,{s})\R^(l,{s})> 



(17) 
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where the projected states \R r (l, {sr}) > and < L r (l, {sl})\ are given by 

w, {«*»> = E4nw.w»> as) 

<^(u^})i = E^n<H r (u^»i (19) 

j ° 

with each state in the summations obtained in a manner analogous to the single-determinantal 
case. The ratio (eqn. (|TTD) is now given as the ratio of sums of determinants appearing 
in the numerator and denominator. Evaluating the ratio hence turns out to be more time 
consuming than in the single determinantal case. 

Property estimates in the single-determinantal PQMC procedure, with a sequential, 
single spin-flip algorithm, are carried out by computing the Green function (eqn. (|15|)) 
at the time slice at which a spin- flip is attempted. This allows the use of an 0(N 2 ) 
updating algorithm for the Green function instead of the usual 0(N 3 ) direct algorithm. 
This is also applicable in the MSPQMC method, except when the states < L J J(l, {sl})| 
and \R^(l, {sr}) > are orthogonal. In this case we use the explicit method of calculating 
matrix elements of the single-particle Green functionpj. The energies presented in this 
communication have been obtained from the Green function (eqn. (fL5|)) estimated at the 
last time-slice, even when a spin flip is attempted at any intermediate time slice. Such 
an estimate of energy will be more accurate, although it precludes the use of an 0(N 2 ) 
updating algorithm. However, we still employ the single sequential spin-flip mechanism 
as it reduces the number of matrix multiplications involved in the computation of the 
Green function. 

In the MSPQMC method for excited states, we encounter the negative sign problem 
even at half-filling although the number of occurrences of the negative signs even at large 
U/t is insignificant (for N = 20, U/t = 6, the fraction of the sample for which negative 
signs are encountered is ~ 10 -4 ). The sign problem here arises because of the phases with 
which the configurations in the trial state are combined although products of individual 
determinants of up and down spin corresponding to < L?J{1, {sL})\Ri r (I, {sr}) > are 
positive. The MSPQMC method, besides being accurate for higher dimensional systems, 
also has the advantage in one-dimensional systems over the DMRG method in that the 
estimates of longer range correlations are as accurate as the nearest neighbour correlations. 

One of the shortcomings of both single- and multi-configurational PQMC calculations 
carried out as described above is that the estimated properties do not reflect the symme- 
tries of the system, as the sampled Ising-configurations do not have the full symmetry of 
the Hamiltonian. To obtain symmetrized property estimates, it is necessary to sample all 
symmetry related Ising- configurations of every Ising configuration that is sampled. We 
have acheived this by a symmetrized PQMC (SPQMC) procedure along the lines of the 
single-determinantal SPQMC met hod 0. Such a symmetrized sampling reduces errors in 
estimates as the sample size is increased by a factor which is of the order of the symmetry 
group, for a marginally small computational overhead. 
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3 Results and Discussion 



The Hubbard chains and rings of < 14 sites can be solved exactly for energies of low-lying 
excitations in various symmetry subspaces. These results provide a strong check on the 
accuracy of the method. Besides, longer Hubbard chains can be solved to a high degree 
of accuracy by the DMRG method and we have compared our MSPQMC results against 
the DMRG results for longer chains (N upto 20). We have computed the energies of the 
lowest excited singlet state connected to the ground state by a dipole transition as well 
the energies of the lowest triplet states for several values of the Hubbard parameter U and 
chain length N. The projection parameter f3 is set to 2.0 with a At of 0.1, all in units 
of t~ x . All the estimates were carried out by averaging over 8000 spin flips per Ising spin 
after allowing 2000 spin flips per Ising spin for equilibration. 

In table 1 we present the MSPQMC energies for three different values of U/t. We note 
that the MSPQMC energies are in very good agreement with exact results (N < 12) or 
high precision DMRG results for (N > 12). The agreement is better at lower correlation 
strengths. While the DMRG energies are better than the MSPQMC energies ||, for the 
chosen states, it is worth noting that the MSPQMC excitation gaps have a better accuracy 
than the absolute energies as the errors in the individual energies for the ground as well 
as the excited states are comparable and have the same sign. 

In fig. (1), we have shown the dependence of the two excitation gaps (with respect 
to the ground state) as a function of U/t for different chain lengths. The 'optical' gap 
increases with U / t while decreasing with N for a fixed U/t. The 'spin' gap decreases with 
U/t with the dependence on iV being similar to the optical gap for a fixed U/t. This 
feature is in agreement with exact as well as DMRG calculations. 

To conclude, we have shown that the PQMC method can be extended to excited states 
using the symmetries of the Hamiltonian via a multi- configurational formulation of the 
PQMC method. Property estimates can be improved by symmetrized sampling along 
the lines of a single-configuration SPQMC procedure. To characterize excited states, we 
have also calculated other quantities such as bond-orders, spin correlations and charge 
correlations. All these quantities are in good agreement with exact/DMRG results and 
will be presented in a longer paper. 

Acknowledgement: We thank Dr. Biswadeb Dutta for help with the computer systems 
at JNCASR and Ms. Y. Anusooya and Mr. Swapan Pati for help with exact and DMRG 
results. 
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Figure Captions 

Figure 1. "Optical" (filled symbols) and "Spin" (open symbols) g function 
of correlation strength (U/t) for Hubbard chains of 16(squares), 18(circles) and 
20 (triangles) sites, from the MSPQMC method. 
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Table 1. MSPQMC energies (A) of lowest dipole allowed singlet and lowest triplet 
excitations (in units of t) for Hubbard chains of length 8 to 20 sites compared with 
exact /DMRG results (B). For chain lengths (> 14 sites) comparison is with 
symmetrized DMRG results with a cut-off of 150. Data for U/t—6.0 is with At = 0.05. 



u/t 


N 


Singlet 


Triplet 






A 


B 


A 


B 


1.0 


6 


-4.4766 


-4.4777 


-4.9331 


-4.9189 




8 


-6.7881 


-6.7906 


-7.1378 


-7.1381 




10 


-9.0185 


-9.0208 


-9.3066 


-9.3083 




12 


-11.2030 


-11.2058 


-11.4470 


-11.4516 




14 


-13.3581 


-13.3634 


-13.5767 


-13.5785 




16 


-15.5007 


-15.5032 


-15.6858 


-15.6948 




18 


-17.6262 


-17.6307 


-17.7970 


-17.8037 




20 


-19.7429 


-19.7494 


-19.8998 


-19.9073 


4.0 


6 


-0.3987 


-0.4221 


-2.7088 


-2.6915 




8 


-1.8982 


-1.9301 


-3.9462 


-3.9165 




10 


-3.2633 


-3.3044 


-5.0796 


-5.1151 




12 


-4.5378 


-4.6062 


-6.2310 


-6.2989 




14 


-5.7782 


-5.8645 


-7.4085 


-7.4734 




16 


-7.0329 


-7.0948 


-8.5557 


-8.6418 




18 


-8.2260 


-8.3059 


-9.6945 


-9.8060 




20 


-9.3868 


-9.5036 


-10.8457 


-10.9672 


6.0 


6 


1.9294 


1.9212 


-1.9793 


-1.9707 




8 


0.7405 


0.6990 


-2.9178 


-2.8677 




10 


-0.3477 


-0.3723 


-3.7362 


-3.7455 




12 


-1.3415 


-1.3639 


-4.6146 


-4.6123 




14 


-2.2349 


-2.3089 


-5.4068 


-5.4724 




16 


-3.1609 


-3.2248 


-6.2614 


-6.3280 




18 


-4.0374 


-4.1214 


-7.1203 


-7.1805 




20 


-4.8692 


-5.0047 


-7.8726 


-8.0308 
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□ 



o 



1 6 sites - optical gap 
16 sites - spin gap 
1 8 sites - optical gap 
18 sites - spin gap 
20 sites - optical gap 
20 sites - spin gap 



i! 



2 







2 



3 



4 



6 



U/t 



